/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */
package testing_class;

import java.util.Random;

/**
 *
 * @author Asko
 */
public class CorrelationFunctions {
    
 public  static  int numOfPoints  ;  //количество реализаций процесса 
 public  static  double arr[] ;
 public  static  Random rand ;
 public static double expected ;
 public static double dispersion ;
 public static double legs[] ;
 public static double t_legs[] ;
 
  
    public static void main(String [] args) {
        
       double gm =  0.879 ;
       double k = 0 , m =  4;
       double omega = 5 ;
       double lambda = 1 ;
       
       legs  = new double[(int)k] ;
       t_legs = new double[(int)k];
      
       int valNum = 5000 ;
       double sampIntrv  =  findSamplingInterval(omega,0.02,2,lambda);
       
        System.out.println("Sample interval -> "  + sampIntrv);
       double th[]  =  theoryFuncN2(13, sampIntrv , lambda);
       double y[]   = realFuncN2(valNum,sampIntrv, lambda);
       double ord[] =  corFuncOrdinates(y, 13,  valNum);
       
    //  double numOrd =  Test.findOrdinatesNum(6,lambda,omega);
  //    double sampIntrv =  CorrelationFunctions.findSamplingInterval(omega, 0.02, 6, lambda);
   
      
    
      
    //    System.out.println("Finding gamma - >  "  + findGamma(5,sampIntrv));
      
      
       //System.out.println("SampInterval - > " + sampIntrv);
     /* double arr[]  = CorrelationFunctions.realFuncN6(valNum, sampIntrv, lambda, omega);
      double ord[]  = CorrelationFunctions.corFuncOrdinates(arr,(int)numOrd , valNum);
     //  System.out.println("---------------->");
      double b[]  = Test.findB(ord,sampIntrv);
        //  System.out.println("---------------->");
      double t[]  = Test.findT(numOrd, sampIntrv);
         // System.out.println("---------------->");
      double a[] =              Test.findA(ord ,b,t );
 
      double Betta[]  =   Test.findBetta(t,a,b,k,gm);
             //System.out.println("-------------------------------------->>>>");
      double res[]  =     Test.approxFuncOrd(Betta,t,gm, ord); */
        
  
        
    } 
    
   
    
    public static double findOrdinatesNum(int modelNum,double lambda, 
             double omega) {
         
        double ordNum =  0 ;  

        switch(modelNum) {

            case 1: { ordNum =  9 ;} break ;      

            case 2: { ordNum = 13 ;} break ;  

            case 3: { ordNum = 19 ;} break ;  

            case 4: { ordNum = (47.0*omega)/(2.0*Math.PI*lambda) ;} break ;  

            case 5: { ordNum = (47.0*omega)/(2.0*Math.PI*lambda); } break ;

            case 6: { ordNum = (47.0*omega)/(2.0*Math.PI*lambda) ;} break ;

            default:   break;    

                            }     

         return Math.round(ordNum) ;
         
     }
    
    public static double findGamma (double k, double sampIntrv) 
    { 
        return 0.4/((2*k + 1) * sampIntrv) ; 
    }
    
    public static double findInverseLegNorm(int k , double gamma) {
        
        
        return  1/2*gamma*(2*k+1);
        
    }
      
    
    
    /*Находим интервал дискретизации*/ 
    public static double findSamplingInterval(double omega ,
            double error, int modelNum,double lambda) {
       
       double derivative2 = 0 ;
       
    switch (modelNum) {
          case 1 : {   derivative2 = lambda * lambda ;  }    break ;    
          case 2 : {  derivative2 =  -lambda * lambda   ; }  break ;               
          case 3 : {  derivative2  = 3 * lambda * lambda ; }  break ;             
          case 4 : {  derivative2  = lambda * lambda - omega*omega ; } break ;                      
          case 5 : {  derivative2  = -lambda * lambda - omega*omega  ; }  break ;         
          case 6 : {  derivative2  =  -lambda * lambda - omega*omega ; }   break ;                  
          default :   break ;    
               
       }
      
   return Math.sqrt(8*error/Math.abs(derivative2));    
        
    }
    
    
    //Моя модель -> e*(cos(w*dT)- l/w*sin(w*dT))
    public static double []  theoryFuncN6(int nY,double sampInterval,
            double lambda,double omega) {
      
      arr  = new double [nY];
      int i = 0  ;
      double add = 0 ;
      
        while(i++<nY) {
            
         arr[i-1] = Math.exp(-lambda*add)*(Math.cos(omega*add)-
                 lambda/omega*Math.sin(omega*add) ) ;  
            
         add  += sampInterval ;   
      //      System.out.println("PY:" + (i-1) +  arr[i-1]);          
        }
       return  arr ;
    }
    
    public static double []  corFuncOrdinates(double [] y, int numOrd,int len)  {
       
      double mean =  0 ,dispers = 0 ; 
      double px[]  =  new double[numOrd];      
      int     i=0 ;
      double sum = 0 ,
             div = 0 ;
        
 /*  while(i++<len) {
       
    System.out.println("Val num["+ i+"] : " + y[i-1]) ; 
          
       
   } */       
  //Считаем матожидание
   while(i<len) {
           
       sum  += y[i]  ;     
       i++  ;  
       } 
   double val = Math.pow(len, -1) ;
       mean =  val * sum ;
       
      // System.out.println("Val value: " + val);
     //  System.out.println("Sum value: " + sum);
       System.out.println("Mean value: " + mean) ; 
       expected = mean ;
       i=0 ;
 //Считаем дисперсию    
   while(i++<len) {
           
       div  += (y[i-1] - mean)*(y[i-1] - mean)  ;     
           
       } 
       dispers = (Math.pow(len-1, -1)) * div ;
       
       dispersion  = dispers ;
       System.out.println("Dispersion: " + dispers +"\n\n") ;    
    
    i = 0;
    
    sum = 0 ;
  for(int j=0;j<numOrd;j++) {
        
    
       for(int ii=0;ii<len-j;ii++)   {
  
         sum += (y[ii]-mean) *(y[ii+j]-mean) ;   
        
        
       }
       
    px[j]   = Math.pow((len-j-1)*dispers,-1) * sum ;  
    sum =  0 ; 
   System.out.println("px[" + j + "] = " + px[j]  ); 
        
    }   
       
        
        
        
        return px ;
    }
      
    public  static double []  realFuncN6(int len ,double sampInterval
            ,double lambda,double omega) {
   
       double x []   = new double[len];
       double y []   = new double[len];
      
       
       rand = new Random();
       
       double alpha0 , alpha1, a0, a1, b1,b2;
       double gamma, gamma0 ;
       double mean, dispers ;
       double sum  = 0 , div = 0;
       
       gamma =  lambda * sampInterval ;
       gamma0 =  omega * sampInterval ;
       
       double p = Math.exp(-gamma);
        
       alpha0 = p*(p*p-1)*Math.cos(gamma0) - lambda/omega*(p*(p*p+1))*Math.sin(gamma0);
       
       alpha1 = 1 - Math.pow(p, 4) + 4*p*p*lambda/omega*Math.sin(gamma0)*Math.cos(gamma0);
       
      
    if (alpha1  >= Math.sqrt( alpha1*alpha1 - 4*alpha0*alpha0  )) {
           
           
       a0 = Math.sqrt((alpha1-Math.sqrt( alpha1*alpha1 - 4*alpha0*alpha0  ))/2) ;   
           
       }  
    else  {
           
       a0 = Math.sqrt((alpha1+ Math.sqrt( alpha1*alpha1 - 4*alpha0*alpha0  ))/2) ;     
           
       }
       
       a1 = alpha0/a0 ; 
       
       
       b1 = 2*p*Math.cos(gamma0);
       b2 = -p*p ;
       
       
       
        int i = 0 ; 
        
    while (i++<len) {
           
     //  System.out.println("Iterator: " +  i )    ;
        x[i-1] =  rand.nextDouble();
   //     System.out.println("Random value number " +  i  + "is"  + x[i-1] );
           
    }
       i = 1 ;
       
           
       y[0] = x[0]*a0;     
       y[1] = a0*x[1] + a1*x[0] + b1*y[0];
       
   while(++i<len) {
           
       y[i] = a0 * x[i] +  a1* x[i-1] + b1*y[i-1] + b2*y[i-2];     
           
   } 
       
 /*     i=0 ;
  //Считаем матожидание
   while(i++<len) {
           
       sum  += y[i-1]  ;     
           
       } 
   double val = Math.pow(len, -1) ;
       mean =  val * sum ;
       
     
       System.out.println("Mean value: " + mean) ; 
       
       i=0 ;
 //Считаем дисперсию    
       while(i++<len) {
           
       div  += (y[i-1] - mean)*(y[i-1] - mean)  ;     
           
       } 
       dispers = (Math.pow(len-1, -1)) * div ;
      
       System.out.println("Dispersion: " + dispers + "\n\n") ; */
    
        return  y ;
    }
    //Третья модель -> e*(1- l/dT)
     public  static double []  theoryFuncN3(int nY,double sampInterval,
             double lambda) {
        
      double [] array  = new double [nY];
      int i = 0  ;
      double add = 0 ;   
        
       while(i++<nY) {
            
         array[i-1] = Math.exp(-lambda*add)*(1- lambda*add ) ;  
            
         add  += sampInterval ;   
        
           System.out.println("T [" + (i-1) +  "] : " +  array[i-1]);   
        } 
        
        
        
        return  array ;
    }
   //Вторая модель -> e*(1 + l/dT) 
      public  static double []  theoryFuncN2(int nY,double sampInterval,
              double lambda) {
        
       double [] array  = new double [nY];
      int i = 0  ;
      double add = 0 ;   
        
       while(i++<nY) {
            
         array[i-1] = Math.exp(-lambda*add)*(1 +  lambda*add ) ;  
            
         add  += sampInterval ;   
        System.out.println("tx [" + (i-1) +  "] : " +  array[i-1]);    
        } 
      return  array ;
    }
     
     
     
     
      public  static double []  realFuncN2( int len, double sampIntrv,
              double lambda) {
        
       double x []   = new double[len];
       double y []   = new double[len];
      
       
       rand = new Random();
       
       double alpha0 , alpha1, a0, a1, b1,b2;
       double gamma, gamma0 ;
       double mean, dispers ;
       double sum  = 0 , div = 0;
       
       gamma =  lambda * sampIntrv ;
      // gamma0 =  omega * sampIntrv ;
       
       double p = Math.exp(-gamma);
        
       alpha0 = Math.pow(p,3)*(1+gamma) - p*(1-gamma);
       
       alpha1 = 1 - (4*Math.pow(p, 2)*gamma)  - Math.pow(p, 4);
       
      
    if (alpha1  >= Math.sqrt( alpha1*alpha1 - 4*alpha0*alpha0  )) {
           
           
       a0 = Math.sqrt((alpha1-Math.sqrt( alpha1*alpha1 - 4*alpha0*alpha0  ))/2.0) ;   
           
       }  
    else  {
           
       a0 = Math.sqrt((alpha1+ Math.sqrt( alpha1*alpha1 - 4*alpha0*alpha0  ))/2.0) ;     
           
       }
       
       a1 = alpha0/a0 ; 
       
       
       b1 = 2*p;
       b2 = -p*p ;
       
       
       
        int i = 0 ; 
        
    while (i++<len) {
           
     //  System.out.println("Iterator: " +  i )    ;
        x[i-1] =  rand.nextDouble();
   //     System.out.println("Random value number " +  i  + "is"  + x[i-1] );
           
    }
       i = 1 ;
       
           
       y[0] = x[0];     
       y[1] = a0*x[1] + a1*x[0] + b1*y[0];
       
   while(++i<len) {
           
       y[i] = a0 * x[i] +  a1* x[i-1] + b1*y[i-1] + b2*y[i-2];     
           
   } 
       
 /*     i=0 ;
  //Считаем матожидание
   while(i++<len) {
           
       sum  += y[i-1]  ;     
           
       } 
   double val = Math.pow(len, -1) ;
       mean =  val * sum ;
       
     
       System.out.println("Mean value: " + mean) ; 
       
       i=0 ;
 //Считаем дисперсию    
       while(i++<len) {
           
       div  += (y[i-1] - mean)*(y[i-1] - mean)  ;     
           
       } 
       dispers = (Math.pow(len-1, -1)) * div ;
      
       System.out.println("Dispersion: " + dispers + "\n\n") ; */
     
        return  y ;
        
    
        
    }
      
      public static double pxFinder(int j , int len) {
          
          
      
          
          
          
          
        return  0 ;  
      }
    
      
//////////////////////////////////////////////////////////////////////////////////      
      
    
     
     public static double leg(double t, double gm , double k) {
         
         
      double res = 0 ; 
      
     if (k == 0)  {
         
      res  = Math.exp(-t*gm);   
         
         
     }   else  { 
         
         for (int i = 0; i < k+1; i++) {
          
  res = res  + factorial(k)/(factorial(k-i)*factorial(i))*
          factorial(k+i)/(factorial(k)*factorial(i)) * 
          Math.pow(-1, i) * Math.exp(-(2*i+1)*gm*t); 
          
          }
     }
        
        
         
      return res ; 
     }
     
     public static double leg0(double t, double gm) {
     
         return Math.exp(-t*gm);
     
     } 
     
     public static double leg1(double t, double gm) {
         
         return Math.exp(-t*gm) - 2*Math.exp(-3* t * gm) ;
         
     }
    //Находим к-ый от времени т
     public static double recurLegFinder(double t, double gm, double k) {
    
      
  if(k == 0)  {
      
          return leg0(t,gm);}
      
   else  
          
          if( k == 1) { 
          return leg1(t,gm);
          }
      
    else  
       return 
             (2*k-1)/k * (1 - 2*Math.exp(-2*gm*t))* recurLegFinder(t,gm,k-1) -
             (k-1)/k * recurLegFinder(t,gm,k-2);

     }
     
     public static double leg0(double t0, double t1, double gm) {
     
         return Math.exp(-t1*gm ) -  Math.exp(-t0 * gm);
     
     } 
   
     
     public static double leg1(double t0, double t1, double gm) {
         
        return  Math.exp(-t1*gm) - 2*Math.exp(-3* t1 * gm) 
               -
                Math.exp(-t0*gm) - 2*Math.exp(-3* t0 * gm)  ;
         
     }
    
     public static double recurLegFinder(double t0, 
             double t1, double gm, double k) {
    
      
  if(k == 0)  {
      
          return leg0(t0,t1,gm);}
      
   else  
          
          if( k == 1) { 
          return leg1(t0,t1,gm);
   
          
          }
      
    else  
       return 
        (2*k-1)/k * (1 - 2*Math.exp(-2*gm*t1))* recurLegFinder(t0,t1,gm,k-1) -
        (k-1)/k * recurLegFinder(t0,t1,gm,k-2)

    -   
        ((2*k-1)/k * (1 - 2*Math.exp(-2*gm*t0))* recurLegFinder(t0,t1,gm,k-1) -
        (k-1)/k * recurLegFinder(t0,t1,gm,k-2))      ;

     }
  
     
     //Интеграл 0 порядка без т  
   public static double legI0(double t, double gm) {
       
       
     return -1/gm*Math.exp(-gm*t);       
       
       
   }  
   //    return  2*Math.exp(-3 * gm * t)/3*gm - Math.exp(- gm * t)/gm ; 
  //Интеграл 1 порядка без т
   public static double legI1(double t , double gm) {
   
     return  -1/gm * (1 - 2.0/3.0 * Math.exp(-2 * gm * t)) * Math.exp(-gm * t );         
       
   } 
   
   //Интеграл 0 порядка с участием т
   public static double legIT0(double t, double gm) {
       
       
     return -(t/gm +  1/(gm*gm))*Math.exp(-gm*t) ;
             
       
       
   } 
  //Интеграл 1 порядка с участием т 
   public static double legIT1(double t, double gm) {
       
       
     return 
             
             -(t/gm +  1/(gm*gm))*Math.exp(-gm*t)
             
             + 
             
             2*(t/(3.0*gm) +  1.0/(9.0*gm*gm))*Math.exp(-3.0* gm * t);
             
           

   } 
  //Рекурсивный поиск интеграла без т
   
   public static  double legRecIntegral(double t , double gm, double k) {
       
       
        double [] arr  =  new double[(int)k+1];  
       
       if (k == 0 )  {  
           
           return legI0(t,gm);
       }
       
       else if (k == 1) {
           
        //   System.out.println("In >>>>");
           return legI1(t,gm);
           
       }
       
       else 
           
           return 2*(2*k-1)/(2*k+1) * legRecIntegral(t,gm,k-1)- 
                   
                  (2*k-3)/(2*k+1) * legRecIntegral(t,gm,k-2) - 
                   
                  1/(gm*(2*k+1))* (leg(t,gm,k) - leg(t,gm,k-2) )  ;
    
   } 
  //Рекурсивный поиск интеграла с участием т
   public static  double legRecIntegralT(double t , double gm, double k) {
       
       
        double [] arr  =  new double[(int)k+1];  
       
       if (k == 0 )  {  
           
           return legIT0(t,gm);
       }
       
       else if (k == 1) {
           
        //   System.out.println("In >>>>");
           return legIT1(t,gm);
           
       }
       
       else 
           
           return  
                   2*(2*k-1)/(2*k+1) * legRecIntegralT(t,gm,k-1) -
                   
                   (2*k-3)/(2*k+1)*legRecIntegralT(t,gm,k-2)  -
                   
                   1/(gm*(2*k+1))  * (t*(leg(t,gm,k) -  leg(t,gm,k-2))  -  
                   
                   legRecIntegral(t,gm,k)+legRecIntegral(t,gm,k-2)   )
                      ;
  
   } 
  //Находим Ai 
   public static double [] findA ( double p[] , double b[] , double t[]) {
       
       double a [] =  new double[b.length];
     
       for (int i = 0; i < a.length; i++) {
       
           a[i]  = p[i] - b[i]*t[i] ;
    //       System.out.println("a[" + i +  "] = " +a[i] );  
       }
       
      return a ;  
   
   }
   //Находим Bi
   public static double [] findB ( double []p, double sampIntrv) {
       
       double b [] = new double [p.length-1];
       
       for (int i = 0; i <p.length-1 ; i++) {
           
          b[i] = (p[i+1] - p[i])/ sampIntrv ; 
   //       System.out.println("b[" + i +  "] = " +b[i] );   
       }

       return b ;  
       
   }
   //Находим массив т-итых
   public static double [] findT (double numOrd , double sampIntrv) {
       
       double t [] = new double[(int)numOrd];
       double iter =  sampIntrv ;
       
       for (int i = 0; i < t.length; i++) {
           
           t[i] = sampIntrv - iter ; 
           
           sampIntrv += iter ; 
     //      System.out.println("T[" + i + "] = "  + t[i]);
       }
      return t;  
       
   }
   //Поиск бетта-коффициентов
   public  static double[] findBetta(double t[] , double a[] ,
           double b[] , double k, double gm) {
    
    
       
       double [] Betta =  new double[(int)(k)];
       double sum = 0 ; 
       
    for(int i=0; i<k; i++) {
        
        for (int j = 0; j < t.length-1; j++) {
            
           sum += a[j]*(legRecIntegral(t[j+1] , gm , i)  -
                   legRecIntegral(t[j] , gm , i) ) 
                   +
                  b[j]*(legRecIntegralT(t[j+1] , gm , i)  - 
                   legRecIntegralT(t[j] , gm , i)); 
                   
                   ; 
            
        }
       Betta[i] = 2*gm*(2*i+1) *  sum; 
       System.out.println("Betta [" +  i  + "]  is - > " + Betta[i] );
       
       sum = 0 ;
        
    }
    
    return Betta ;
    
   }
   
  
   
   public static double [] approxFuncOrd(double Betta[] , double ti[]
           , double gm, double px[]) {
       
       
      double []  Kx = new double[ti.length];
      double sum =  0 ;
  
  for (int i = 0; i < ti.length; i++) {
           
       
   
      for(double j = 0 ; j<Betta.length;j++) {
          
          
      sum +=Betta[(int)j]*leg(ti[i],gm,j);  
          
          
           }
      Kx[i] = sum ;
    System.out.println("Kx[" + i + "] = " + Kx[i]  + "   " +
            " px[" + i + "] = " +  px[i]);   
              
              sum = 0 ;
              
       }   
  
      return Kx ; 
       
       
   }
   
   public static double [] findPowerSpectrum(double Betta[], double gm 
           , double omega, double step) {
      
       int length =    (int) (2*(omega / step)); 
       
       
      double iter =  omega ; 
      double S[]   =   new double[length];
      double sum   =  0 ;  
      double FiSum =  0 ;
      omega = 0 ;
      
        for (int i = 0; i < S.length;i++) {
            
           
            
                  for (int j = 0; j < Betta.length; j++) {
                      
                      
                      if (j==0)  {
                          
                         FiSum = 0 ; 
                      }
                      
                      else {
                          
                          for (int l = 0; l < j; l++) {
                              
                           FiSum += Math.atan(omega/(2*gm*(l+1))) ;   
                              
                          }
             
                      }
                      
                   sum += Betta[j] *   (1.0/(2*j+1)) * 
                       Math.cos(Math.atan(omega/(2*gm*(j+1))))  * 
                       Math.cos(Math.atan(omega/(2*gm*(j+1))) + 2* FiSum  ) ;
                       
                   FiSum = 0 ;    
                      
                      
            }
                  
          S[i]  =  (1.0/(Math.PI*gm))*sum ;       
            System.out.println("S [" + i + "]  - > "  +  S[i]);     
            
          
          sum  = 0 ;
          omega += step ;
          
           

        }
    
      return S ; 
       
   }
     
    
     public static double []  findTheorySpectrumN2(double lambda, double omega,
             double step) {
       
         int len =   (int)(2*(omega/step));
         double om = 0 ;
         double [] s = new double [len];
         
         for (int i = 0; i < s.length; i++) {
 
             s[i]  =  (2.0 * Math.pow(lambda, 3)) / 
                     (Math.PI*(Math.pow((lambda*lambda +  om*om), 2))); 
             
             om += step ;
             
         }
         
         
         return s ;
         
         
     }
    
     
     public static double []  findTheorySpectrumN6(double lambda, double omega,
             double step) {
       
         int len =   (int)(2*(omega/step));
         double om = 0 ;
         double [] s = new double [len];
         
         for (int i = 0; i < s.length; i++) {
 
             s[i]  =  (2.0*lambda*om*om) / 
                     (Math.PI*
                   
                     (Math.pow(lambda, 2) + Math.pow(om - omega,2)) *
                     
                     (Math.pow(lambda, 2) + Math.pow(om + omega,2))) ; 
             
             om += step ;
             
         }
         
         
         return s ;
         
         
     }
     
     public static double  factorial(double val) { 
     if (val == 0) return 1 ;
     else  return val*factorial(val-1);
     
     }  
    
}
